#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>
#define k 0.001
#define N 1000000

int func (double t, const double y[], double f[], void *params)
{
	(void)(t); 
	double mu = *(double *)params;
	f[0] = y[1];
	f[1] = -(2/t) * y[1]-mu*mu*y[0];
	return GSL_SUCCESS;
}

int	jac (double t, const double y[], double *dfdy, double dfdt[], void *params)
{
	(void)(t);
	double mu = *(double *)params;
	gsl_matrix_view dfdy_mat = gsl_matrix_view_array (dfdy, 2, 2);
	gsl_matrix * m = &dfdy_mat.matrix;
	gsl_matrix_set (m, 0, 0, 0.0);
	gsl_matrix_set (m, 0, 1, 1.0);
	gsl_matrix_set (m, 1, 0, -2.0*mu*mu*y[1]/t);
	gsl_matrix_set (m, 1, 1, -2.0*mu*mu*y[0]/t);
	dfdt[0] = 0.0;
	dfdt[1] = 2*y[1]/(t*t);
	return GSL_SUCCESS;
}

int main (void)
{
	double mu = k;
	gsl_odeiv2_system sys = {func, jac, 2, &mu};
  
	gsl_odeiv2_driver * d = gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk8pd, 1e-6, 1e-6, 0.0);

	double t = 1e-2, t1 = 1e4;
	double y[2] = { 1.0, 0.0 };

        char *path = "result.txt";
        FILE* filepath = fopen(path, "w");
	for (int i = 1; i <= N; i++)
	{
    	double ti = t + (i * (t1 - t) / N);
    	int status = gsl_odeiv2_driver_apply (d, &t, ti, y);

    	if (status != GSL_SUCCESS)
        {
        	printf ("error, return value=%d\n", status);
        	break;
        }

//    	printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
    	fprintf(filepath,"%.5e %.5e %.5e\n", t, y[0], y[1]);
    }

	gsl_odeiv2_driver_free (d);
	return 0;
}
